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ABSTRACT 

The initial conditions and relevant physics for the formation of the earliest galaxies are well specified 
in the concordance cosmology. Using ab initio cosmological Eulerian adaptive mesh refinement radia- 
tion hydrodynamical calculations, we discuss how very massive stars start the process of cosmological 
reionization. The models include non-equilibrium primordial gas chemistry and cooling processes 
and accurate radiation transport in the Case B approximation using adaptively ray traced photon 
packages, retaining the time derivative in the transport equation. Supernova feedback is modeled by 
thermal explosions triggered at parsec scales. All calculations resolve the local Jeans length by at 
least 16 grid cells at all times and as such cover a spatial dynamic range of ~10^. These first sources 
of reionization are highly intermittent and anisotropic and first photoionize the small scales voids 
surrounding the halos they form in, rather than the dense filaments they are embedded in. As the 
merging objects form larger, dwarf sized galaxies, the escape fraction of UV radiation decreases and 
the H II regions only break out on some sides of the galaxies making them even more anisotropic. In 
three cases, SN blast waves induce star formation in overdense regions that were formed earlier from 
ionization front instabilities. These stars form tens of parsecs away from the center of their parent DM 
halo. Approximately 5 ionizing photons are needed per sustained ionization when star formation in 
10^ Mq halos are dominant in the calculation. As the halos become larger than ~lO''Af0, the ionizing 
photon escape fraction decreases, which in turn increases the number of photons per ionization to 
15-50, in calculations with stellar feedback only. Radiative feedback decreases clumping factors by 25 
per cent when compared to simulations without star formation and increases the average temperature 
of ionized gas to values between 3,000 and 10,000 K. 

Subject headings: cosmology: theory — intergalactic medium — galaxies: formation — stars: forma- 
tion 



1. MOTIVATION 

It is clear that quasars are not responsible to keep 
the universe ionized at redshift 6. The very brightest 
galaxies at those redshifts alone also provide few pho- 
tons. The dominant sources of reionization so far are 
observationally unknown despite remar kable advances 
in finding sources at high redshift (e .g. 'Shapircj Il986j: 
Bouwens et all 120041: iFan et al.l l200d [Thompson et al l 



20071: lEvles et al.ll2006l) and hints for a la rge number of 



unres olved sources at very high redshifts (jSpergel et al.l 
120071 : ^Kashlinsky et al. 200 3) which is still a top ic of 
debate (|Coorav et al.l 120071 : iThompson et al.ll2007h . At 
the same time, ab initio numerical simulations of struc- 
ture formation in the concordance model of structure 
formation have found that the first luminous objects 
in the universe are formed inside of cold dark matter 
(CDM) dominated ha l os of total mass e s 2 x 1 0^ - lO^MfD 
(iHaiman et all 119961 : iTegmark et"an 119971 : lAbel et all 
1998f). Fully cosmol ogical ab initio c alculations of 



Abel et al.l (|2000ll2002[ ) and more recently lYoshida et all 
( 20061 ) ciearlv show that these objects will form isolated 



very massive stars. Such stars will be copious emitters of 
ultraviolet (UV) radiation and are as such prime suspect 
to get the process of cosmological reioni zation started. 
In fac t, on e dimensional c a lculat ions of IWhalen et al.l 
(|2004( ) and iKitavama et al.l (|2004f ) have already argued 
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that the earliest H II regions will evaporate the gas from 
the host halos and that in fact most of the UV radiation 
of such st ars would escape into t he int ergalactic in e dium 
Recently, lYoshida et all (j2007aD and lAbel et all (|2007f ) 
demonstrated with full three-dimensional radiation hy- 
drodynamical simulations that indeed the first H II re- 
gions break out of their host halos quickly and fully 
disrupt the gaseous component of the cosmological par- 
ent halo. All of this gas finds itself radially moving 
away from the star at ~ 30 km s~^ at a distance of 
~ 100 pc at the end of the stars life. At this time, 
the photo-ionized regions have now high electron frac- 
tions and little destructive Lyman- Werner band radia- 
tion fields creating ideal conditions for molecular hydro- 
gen formation which may in fact stimulate further star 
formation above levels that would have occurred without 
the pre-ionization. Such conclusion have been obtained 
in calculations with approximations to multi dimen- 
sional r adiative tra nsfer or one dimensional numerical 
models (iRicotti eTal. 2002a; Nagakura fc Om ukai 200^; 
O'Shea et al.ri2 005: Yos hida et al.ll200'6l : lAhn fc Shapirol 



20071 : iJohnson et al.. ,20071 ). These early stars may also 



explode in supernovae and rapidly enrich the surround- 
ing material with heavy elements, deposit kinetic energy 
and entropy to the gas out of which subsequent struc- 
ture is to form. This illustrates some of the complex 
interplay of star formation, primordial gas chemistry, ra- 
diative and supernova feedback and readily explains why 
any reliable results will only be obtained using full ab 
initio three dimensional hydrodynamical simulations. In 
this paper, we present the most detailed such calculations 
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yet carried out to date and discuss issues important to 
the understanding of the process of cosmological reion- 
ization. 

It is timely to develop direct numerical models of early 
structure formation and cosmological reionization as con- 
siderable efforts are underway to 

1. Observationally find the earliest galaxies with 
the James Webb Space Telescope (JWST; 
[Gardner et al.l 12006^ an d the At acama Large 
Millimeter Array (ALMA: iWilson e t al. 2005i), 

2. Further constrain the amount and spatial non- 
uniformity of the polarization of the cosmic m i- 
crowave background radiation (jPage et al.ll2007l ). 

3. Measure the surface of reionizatioi i with LOFAR 
(Rottgering et al. 20 Q1), MWA (jBowman et all 
[2 007 )~GMRT (Swarup et al/'iggf) and the Square 
Kilometer Array (SKA: iSchilizzi 2004 ). and 

4. Find high redshift g amma ray bursts with SWIFT 
(jOehrels et alJ[20M ) and their infrared follow up 
observations. 

We begin by describing the cosmological simulations 
that include primordial star formation and accurate ra- 
diative transfer. In we report the details of the star 
formation environments and host halos in our calcula- 
tions. Then in 21 "^6 describe the resulting start of cos- 
mological reionization, and investigate the environments 
in which these primordial stars form and the evolution 
of the clumping factor. We compare our results to pre- 
vious calculations and further describe the nature of the 
primordial star formation and feedback in ^ Finally we 
summarize our results in the last section. 

2. RADIATION HYDRODYNAMICAL 
SIMULATIONS 

We use radiation hydrodynamical simulations with 
a modified version of the cosmological AMR code 
Enzo to study the radiative effec ts from the first 
stars (jBrvan fc NormanI Il997l 1X9991). We have inte- 
grated adaptive ray tracing ( jAbel fc Wandeltll2002 f) into 
the chemistry, energy, and hydrodynamics solvers in 
Enzo that accurately follow the evolution of the H II 
regions from stellar sources and their relevance during 
structure formation and cosmic reionization. 

Seven different simulations are discussed here. Table [T] 
gives an overview of the parameters and the physics in- 
cluded in these calculations. We perform two cosmo- 
logical realizations, Sim A and B, with three sets of as- 
sumptions about the primordial gas chemistry. The sim- 
plest calculations here assume only adiabatic gas physics 
and provide the benchmark against which the more in- 
volved calculations are compared. We compare this to 
one model with atomic hydrogen and helium cooling only 
and one that includes H2 cooling. Massive, metal-free 
star formation is included only in the II2 cooling models. 

These calculations are initialized at redshift'^ z = 130 
(120) when the intergalactic medium has a temperature 
of 325 (280) K in box sizes 1 comoving Mpc (1.5 Mpc) for 

^ To simplify the discussion, simulation A will always be quoted 
first with the value from simulation B in parentheses. 



Sim A (B) with three initial nested AMR grids. We use 
the COSMICS package to create the initi al conditions, 
which uses a Zel'dovich approximation (jBertschingeil 
[19951 [200lh . We use the cosmological parameters of 
(^2iJ/l^ Hm, h, as, n) = (0.024, 0.27, 0.72, 0.9, 1) 
from first year WM AP results, where the constants have 
the usual meaning (jSpergel et "al1 l2003f). The ch anges 
in the third year WMAP results (jSpergel et all [20OT ) 
does not affect the evolution of individual halos stud- 
ied here but only de lays structure formation by ~40% 
varez et aLll2006bl ). The adiabatic simulations as well 
as the atomic hydrogen and helium cooling only calcu- 
lations are described in I Wise fc Abel (|2007a^ . The new 
models presented here have the exact same setup and 
random phases in the initial density perturbation and 
only differ in that they include star formation as well 
as follow the full radiation hydrodynamical evolution of 
the H II regions and supernova feedback in Sim B. We 
use the designations RT and SN to distinguish cases in 
which only star formation and radiation transport were 
included (RT) and the one model which also includes su- 
pernovae (SN) in Sim B. We use the same refinement 
criteria as in our previous work, where we refine if the 
DM (gas) density becomes three times greater than the 
mean DM (gas) density times a factor of 2', where / is 
the AMR refinement level. We also refine to resolve the 
local Jeans length by at least 16 cells. Cells are refined to 
a maximum AMR level of 12 that translates to a spatial 
resolution of 1.9 (2.9) comoving parsecs. This spatial 
resolution of ~0.1 proper pc is required to model the 
formation of the D-type front at small scales correctly. 
Refinement is restricted to the innermost initial nested 
grid that has a side length of 250 (300) comoving kpc. 

These simulations focus on a highly biased region at 
z > 15, which could be subject to numerical artifacts 
created by incorrect second- and higher-order growing 
modes (i.e. tra nsients) associated with the Zel ' dovic h 
approximation ([Scoccimarr ol fl998l : iCrocce et al.l[2006h . 
This suppresses the tails of the density and velocity prob- 
ability distributions that leads to rare events (i.e. ha- 
los) becoming less common than expected and is usu- 
ally avoided by starting at very high redshifts. Although 
this may be the case, the halos examined in our sim- 
ulations are well within cosmic variance. To illustrate 
this, we plot the DM halo mass function of the entire 
simulation volume and compare it to the one computed 
by a n ellipsoidal variant (S-T) of Press-Schechter for mal- 
ism (|Press &: Schechtedll974HSheth fc TormenI 12002) in 
Figure [TJ The circles represent the data and its associ- 
ated Poisson errors, and the lines are the S-T number 
densities. The halo mass function matches well with the 
S-T function (solid fine) above ~ lO^-^Mg. Below this 
mass, it decreases relative to the PS function because 
our high-resolution region samples a central box with a 
side L = 250 (300) comoving kpc. To correct for this 
effect, we multiply the S-T mass function by {L/1 Mpc)^ 
to obtain the dotted line. We chose this region because 
it is highly biased and should expect the halo mass func- 
tion to be greater than the cosmic average. The halo 
mass function below ~ lO^'^M© matches well with S-T 
function corrected by a factor of b{L/l Mpc)"^, where b 
= 35 (20) is the overdensity of halos in the refined re- 
gion, that is plotted as the dotted line in Figure [TJ We 
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TABLE 1 

Simulation Parameters 



Name 


I 

[Mpc] 


Cooling model 


SF 


SNe 


Npart 




Nccll 


SimA-Adb 


1.0 


Adiabatic 


No 


No 


2.22 X lO'' 


30230 


9.31 X 10^ (453^) 


SimA-HHe 


1.0 


H, He 


No 


No 


2.22 X lO'' 


40601 


1.20 X 10** (4943 ) 


SimA-RT 


1.0 


H, He, H2 


Yes 


No 


2.22 X 10^ 


44664 


1.19 X 10** (4933) 


SimB-Adb 


1.5 


Adiabatic 


No 


No 


1.26 X 10^ 


23227 


6.47 X 10'^ (4023) 


SimB-HHe 


1.5 


H, He 


No 


No 


1.26 X 10'' 


21409 


6.51 X lO'^ (4023) 


SimB-RT 


1.5 


H, He, H2 


Yes 


No 


1.26 X lO'' 


24013 


6.54 X lO'' (403^) 


SimB-SN 


1.5 


H, He, H2 


Yes 


Yes 


1.26 X lO'' 


24996 


6.39 X lO'^ (400^) 



Note. — Col. (1): Simulation name. Col. (2): Box size. Col. (3): Cooling model. Col. (4): Star formation. Col. (5): Supernova 
feedback. Col. (6): Number of dark matter particles. Col. (7): Number of AMR grids. Col. (8): Number of unique grid cells. 



note that recent high-resolution simulations have found 
that the S-T halo mass function overestimates rare ob- 
jects by up to 50% at all times and halo abundances 
at very high redsh ifts (|Iliev et al.ll2006l: iReed et al.ll2007t 
iLukic et al.ll2007[ ). Transients from the Zel'dovich ap- 
proximation used in our initial setup should not affect 
our general conclusions on the star formation rate and re- 
sulting reionization; however we caution that the simula- 
tions presented here probably underestimates the clump- 
ing factor in the intergalactic medium (IGM) and number 
densities of low-mass halos. The effects of transients in 
simulations that focus highly biased regions should be 
examined more carefully in the future. 

The star formation rec i pe and radiation transport are 
detailed in lWise fc Abell |2007c[ ). Here we overview the 
basics about our method. Star formation is modelled us- 
ing the ICen &: Ostriked (|1992 ') algorithm with the addi- 
tional requirement that an H2 fraction of 5 x 10^'* must 
exist before a star forms. We allow star formation to 
occur in the Lagrangian volume of the surrounding re- 
gion out to three virial radii from the most massive halo 
at 2: = 10 in the dark matter only runs as discussed 




Fig. 1. — Halo mass functions for Simulation A at 2; = 15.9 (left) 
and B at 2: = 16.8 (right). The open circles are computed from 
the simulations. The solid line is the halo mass function from the 
ellipsoidal variant of Press-Schechter formalism l lSheth fc TormenI 
[2Mal . The dotted line is corrected by a factor of (L/1 Mpc)^, 
where L = 250 (300) kpc is the comoving length of the inner refined 
region. The dashed line is corrected by a factor of b(L/l Mpc)^, 
where fe = 35 (20) is the overdensity of halos in the biased, refined 
region where we allow stars to form. Our simulations agree well 
with the uncorrected S-T function (solid) above ~ 10^'^Mq and 
the corrected function (dashed) below this mass scale. 



in IWise &: Abe] (|2007al ). This volume that has a side 
length of 195 (225) comoving kpc at z = 30 and 145 
(160) comoving kpc at the end of the calculation. The 
calculations with SNe use a stellar mass of 170Mq, 
whereas the ones without SNe use a mass of IOOM0. The 
ionizi ng luminosities are taken from no mass loss mod- 
els of ' Schaerer' (! 200^. and w e employ the SN energies 
from lHeger fc WooslevI ()2002D . Star particles after main 
sequence are tracked but are inert. There is evidence 
of lower mass prim ordial stars forming within relic H II 
regions ( jO'Shea et al. 2005; Yoshida et al. 2007b), but 
we neglect this to avoid additional uncertain parameters. 
This is a desired future impro vement, however. 

We use adaptive ray tracing ()Abel fc Wandeltj[2002f) to 
calculate the photo-ionization and heating rates caused 
by stellar radiation. We consider photo-ionization from 
photons with an energy of 28.4 (29.2) eV that is the 
mean energy of ionizing radiation from a metal-free star 
with 100 (170) Mq. The simulation box is large enough 
so that the H II region never expands outside the box; 
however, our ray tracing code employs isolated boundary 
conditions so that photons are lost if they travel outside 
the computational domain. We account for H2 photo- 
dissociation with a 1/r^ Lyman- Werner radiation field 
without self-shielding. We use a non-equilibrium, nine- 
species (H, H+, He, H e+, He++, e" H2, H+, H") 
chemistry solver in Enzo (jAbel et al.ll997HAnninos et al.l 
Il997f ) that takes into account the additional spatial de- 
pendence of the photoionization rates provided by the 
radiation transport. 

We end the simulations when the most massive halo 
begins to rapidly collapse (i.e. icooi < ^dyn) in the hydro- 
gen and helium cooling only runs at redshift 15.9 (16.8). 
The virial temperature Tvir of the halo is ^10"* K at these 
redshifts. 

3. STAR FORMATION 

Here we describe the aspects of massive metal-free 
star formation in our simulations. The first star forms 
at redshift 29.7 (30.8) in halo typical of Pop HI star 
formation witho ut any feedback that has a mass of 
^ 5x lO^MfD (cf.lAb cl et al. 2000|, [200I iMadracekejalJ 
I2OOII : lYoshid a et al. 2003, 2006). Afterwards there are 
a total of 19, 29, and 24 instances of star formation in 
SimA-RT, SimB-RT, and SimB-SN, respectively. 

3.1. Star Formation Rate 

We show the cumulative star formation rate (SFR) in 
units of comoving Mq Mpc"'^ in Figure [2l This quan- 
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TABLE 2 

Selected Star Forming Halo Parameters 



# 


Sim 


Type 


z 


Mvir 




Pc 












L W J 




[cm^'^l 

L J 


L J 


1 


SimB-RT 


1 


30.9 


4.7 X 10^ 


0.081 


1600 


340 


2 


SimA-RT 


1 


29.9 


4.8 X 10^ 


0.094 


6500 


350 


3 


SimB-RT 


4 


23.7 


5.3 X 10^ 


0.059 


2400 


410 


4 


SimB-SN 


4 


21.0 


1.1 X lO'^ 


0.045 


1800 


480 


5 


SimA-RT 


4 


20.4 


6.3 X 10** 


0.069 


550 


440 


6 


SimB-SN 


2 


20.1 


2.6 X 10** 


0.12 


120 


390 


7 


SimB-RT 


2 


19.9 


2.8 X 10*5 


0.12 


870 


450 


8 


SimA-RT 


2 


19.3 


2.9 X 10*5 


0.13 


1300 


440 


9 


SimB-SN 


3 


19.3 


2.3 X 10^ 


0.12 


360 


450 


10 


SimB-RT 


5 


16.8 


3.1 X 10^ 


0.089 


4100 


2500 


11 


SimB-SN 


5 


16.8 


2.9 X 10^ 


0.065 


1100 


590 


12 


SimA-RT 


5 


16.1 


3.0 X 10^ 


0.061 


130 


470 



Note. — Col. (1): Halo number. Col. (2): Simulation source. Col. (3): Star formation type. Col. (4): Redshift. Col. (5): Virial mass. 
Col. (6): Baryon mass fraction. Col. (7): Central number density. Col. (8): Central temperature. 



tity is simply calculated by taking the total mass of stars 
formed at a given redshift divided by the comoving vol- 
ume where stars are allowed to form (see 311) ■ In this 
figure, we decrease the SFR of SimB-SN by a factor of 
1.7 in order to directly compare the rates from the other 
two simulations. This minimizes some of the uncertain- 
ties entered into our calculations when we chose the free 
parameter of Pop III stellar mass. The cumulative rates 
are very similar in both realizations. The refined volume 
of Sim A (Sim B) has an average overdensity 5 = p/ p = 
1.4 (1.8). The more biased regions in Sim B allows for 
a higher density of star-forming halos that leads to the 
increased cumulative SFR. 

We also overplot the cu mulative SFR in atomic hy - 
drogen cooling halos from iHernquist fc Springel ()2003l ) 
in this figure. It is up to an order of magnitude lower 
than the rates seen in our calculations up to redshift 20. 
They only focused on larger mass halos in their simula- 
tions. The disparity between the rates is caused by our 
simulations only sampling a highly biased region, where 



we focus on a region containing a 3-cr density fluctua- 
tion, a nd from the contribution from Pop III stars. The 
rates of lHernauist fc Sprineell pOOSh are calculated from 
an extensive suite of smoothed particle hydrodynamics 
simulations that encompasses both large and small sim- 
ulation volumes and give a more representative global 
SFR due to their larger sampled volumes. However, our 
adaptive spatial resolution allows us to study both the 
small- and large-scale radiative feedback from Pop III 
stars, which is the main focus of the paper, in addition 

to the quantitative measures such as a SFR; 

To estimate a volume-averag ed SFR (i.e. iMadau et all 
Il996l in units of comoving Mq yr~^ Mpc~'^) from the cu- 
mulative SFR, we need to smooth the discretely increas- 
ing cumulative SFR to ensure its time derivative (i.e. 
the SFR) is a smoothly varying function. Otherwise, the 
SFR would be composed of delta functions when each 
star forms. We first fit the cumulative SFR with a cu- 
bic spline with 10 times the temporal resolution. Then 
we smooth the data back to its original time resolution 
and evaluate its time derivative to obtain the SFR that 
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Fig. 2. — Cumulative star formation rate in units of comoving 
Mq Mpc-3 of SimA-RT (solid), SimB-RT (dashed), and SimB-SN 
(dot-dashed). The star formation rate of SimB-SN has been scaled 
by 100/170, which is the ratio of Pop IH stellar masses used in 
SimB-RT and SimB-SN, in order to make a direct comparison be- 
tween the two simulations. The dot-dot-dashed line represents the 
cum ulative s tar formation rate in atomic hydrogen cooling halos 
from IHernquist Springet (i2003i) . 
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Fig. 3. — Comoving star formation rate in units of Mq yr ^ 
Mpc~^. The lines representing the simulation data have the same 
meaning as in Figure |2] The crosses show at which redshifts stars 
form. The rates in SimB-SN are scaled for the same reason as in 
Figure [2] For comparison, we overplot the star formation rates 
from Hernquist & Springel (2003) in atomic hydrogen line cooling 
halos and .Yoshida et al.. (,2003i) for lOOM© Pop III stars. 
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Fig. 4. — Star formation times versus host halo DM mass for 
SimA-RT (circles), SimB-RT (squares), and SimB-SN (triangles). 
One symbol represents one star. The numbers correspond to the 
halo numbers listed in Table [2] 



we show in Figure O We also mark the redshifts of star 
formation with c rosses. We agai n compare our rates to 
ones calculated in lHernguist fc Sp ringcl (2003) for metal- 
enriched stars and lYoshida et aP (|2003i) for Pop III stars 
with a mass of 100 Mq. Our rates are higher for reasons 
discussed previously. We do not advocate these SFRs as 
cosmic averages but give them as a useful diagnostic of 
the performed simulations. 

We see an increasing function from 5 x 10^"' Mq yr^^ 
Mpc-3 at redshift 30 to -6 x lO'^ Mq yr-^ Mpc'^ at 
redshift 20. Here only one star per halo forms in ob- 
jects with masses <, 5 x 10® Mq. Above this mass scale, 
star formation is no longer isolated in nature and can 
be seen by the bursting nature of the star formation af- 
ter redshift 20, wher e the SFR fluctuate s around 10^^ 
Mq yr-i Mpc-3 (cf. iRicotti et"aI]l2002bD . Since we ne- 
glect H2 self-shielding, the strong Lyman- Werner (LW) 
radiation dissociates almost all H2 in the host halo and 
surrounding regions. Thus we rarely see simultaneous in- 
stances of star formation. However, the regions that were 
beginning to collapse when a nearby star ignites form a 
star 3-10 million years after the nearby star dies. This 
only results in a minor change in the timing of star for- 
mation. Furthermore this delay is minimal compared to 
the Hubble time and does not affect SFRs. 

3.2. Star Forming Halo Masses 

We show the star formation times versus the host halo 
DM masses as a function of redshift in Figure |4l The 
DM halo masses are calc ulated with the HOP algorithm 
pisenstein fc Hutlll998f ). First we focus on star forma- 
tion in the largest halo. Around redshift 30, the first star 
forms in all three simulations with a mass of ^5 x IO^Mq . 
The stellar radiation drives a ~30 km s^^ shock wave 
that removes almost all of the gas from the shallow poten- 
tial well. It takes approximately 75 (40) million years for 
gas to reincorporate into the potential well from smooth 
IGM accretion and mergers. At z ^ 24 in SimB-RT, 
the second star forms in the most massive progenitor 
that now has a mass of 4 x IO^Mq. In SimA-RT, the 
merger history is calmer at z = 24 — 30, and enough gas 
becomes available for H2 cooling and star formation at 
z ^ 20. Here the second star forms in the most massive 



progenitor that has a mass of 5 x 10®A/q. In both RT 
simulations, the stellar feedback expels most of the gas 
from its host once again. For Sim A (Sim B), another 
10 (30) million years passes before the next star to form 
in this halo. Once the halo has a mass of ~10^ Mq, the 
potential energy is great enough to confine most of the 
stehar and SNe outflows. In SimA-RT and SimB-RT, 
halos above this mass scale host multiple sites of star 
formation that is seen in the nearly continuous bursts 
of star formation in the most massive halo. SimA-RT 
forms stars more intermittently than SimB-RT because 
it under goes two major merg ers between redshift 17 and 
21 (see iWise fc Abdl l2007a[) . In SimB-SN at z = 21, 
three stars form in succession in the most massive halo. 
Their aggregate stellar and SNe feedback expels the gas 
from its halo one more time. This halo only forms an- 
other star at z = 16.9 (55 million years later) in this halo 
when enough gas has been reincorporated.** 

Most of the stars form in low-mass halos with masses 
^10^ Mq that are forming its first star between redshifts 
18 - 25 in our calculations. A slight increase in host 
halo masses with respect to redshift mainly occurs be- 
cause of the negative feedback from photo-evaporation 
of low -mass halos that are close to other star-forming 
halos (jHaiman et al.l [20011 ). Additional delays in star 
formation may be caused by ultravio let heating and 
H2 dissociation front previous stars (e.g.lMachacek et al.l 
1200 It lYoshida et all l2003t iMesigner et al.l |2006[) . which 
increase the critical halo mass in which gas can cool and 
condense. 

One interesting difference in SimB-SN from the other 
calculations is that star formation is sometimes induced 
in overdensities within the same halo when a SN blast 
wave overtakes it. This occurs in three halos with masses 
of 2.0, 1.5, and 3.8 xIO^Mq at redshifts 19.9, 19.1, and 
17.8, respectively. The same halos in SimB-RT do not 
form two stars before their gas are expelled and thus 
quenching subsequent star formation. 

3.3. Star Formation Environments 

We further study the nature of high-redshift star 
formation by selecting four star forming regions from 
each simulation and studying the surrounding interstellar 
medium (ISM) prior to star formation. T he ISM in the 
10^ K h alos are described in more detail in lWise fc Abell 
(|2007cf ). The sample of regions are chosen in order to 
compare different star formation environments. These 
regions can be categorized into (1) first star inside an 
undisturbed halo, (2) first star that is delayed by LW ra- 
diation, (3) induced star formation by positive feedback, 
(4) star formation after gas reincorporation, and (5) star 
formation in a halo with a virial temperature over 10"* K. 
The represented halos and their parameters are listed in 
Table [2] and annotated in Figure HI 

We plot the mass-weighted radial profiles of number 
density (left columns) and temperature (right columns) 
within the virial radius for these twelve halos in Figure 
[5] and describe them below. 

1. First star (Halo 1, 2) — These stars are the first 
to form in their respective simulation volume. The 

We have run SimB-SN past z = 16.8 and have seen that it 
starts to host multiple sites of star formation. 
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Fig. 5. — Radial profiles of number density (left column) and temperature (right column) for selected star forming halos inside the virial 
radius. We ovcrplot the radially averaged DM density (solid line) in units of cm^"^. These data represent the state of the region 
immediately after star formation. Notice the added complexity (range) in the density and temperature with increasing host halo mass, 
especially if the region has been affected by stellar radiation, as in Halos 3, 4, 5, 10, 11, and 12. The occasional discontinuities in density 
in the inner parsec arise from our star formation recipe when we remove half the mass in a sphere containing twice the stellar mass. 
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Fig. 6. — Density-squared weighted projections of gas density 
(left) and temperature (right) of Sim A. The field of view is 8.5 
proper kpc (1/216 of the simulation volume) and the color scale is 
the same for all simulations. 



structure of the host halos within our resolution hmit 
exhibit similar characteristics, e.g., a self-similar col- 
lapse and c entral temperatu r es of 300 K, as in previ- 
ous studies (lAbel et al'l 12001 l2002t iBromm et aLlbooa 
lYoshidaet al.ll2006[ ). The halo masses are 4.8(4.7) x 10^ 
M(7). Heating from virialization raises gas temperatures 
to 3000 K, and in the central parsec, H2 cooling becomes 
effective and cools the gas down to 200 K that drives the 
further collapse. The mass-weighted central gas densities 
and temperatures are approximately 3000 cm^'^ and 320 
K, respectively. 

2. Delayed first star (Halo 6, 7, 8) — The host ha- 
los have similar radial profiles as the halos that hosted 
the first stars but with masses of >,1O^M0. Here the 
H2 cooling has been stifled by the LW radiation from 
nearby star formation. Only when the halo mass passes 
a critical mass, the core can cool and conden se by 
H2 formation (M achacc k et al.ll200lUYoshida et al.ll2003l : 
lO'Shea fc Nomianii2007nWise fc Abe]|l2007bn . The cen- 
tral densities are lower than the first stars with 1300, 870, 
and 120 cm-3 in SimA-RT, SimB-RT, and SimB-SN, 
respectively. The central temperatures are marginally 
higher at 440, 450, and 390 K. 





3. Induced star formation (Halo 9) — At 



19.3, 



Fig. 7. — Same as Figure [6] but for Sim B. 



a massive star explodes in a SN, whose shell initially 
propagates outward at 4000 km s^^. After 160 kyr, 
the shell passes an overdensity within the same halo 
that is caused by an ioniza tion front instability (e.g., see 
IWhalen fc Norman|[2008allb[ ) . The star forms 35 pc away 
from the SN explosion and the DM halo center. The com- 
bination of the shock passage and excess free electrons in 
the r elic H II cataly z e H2 formation in th i s low-mass halo 
(e£. iFerraral [19981: lO'Shea et all 120051 : iMesigner etall 
I2006f) . The SN blast wave heats the gas over 10^ K to 
radii as 1 pc. In the density profile, both low and high 
density gas exists at similar radii. Here the shock pas- 
sage creates a tail of gas streaming from the central core, 
whose asymmetries can be seen in the density profile. 
However, the core survives and benefits from the excess 
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electrons created during this event. The central tem- 
perature is similar to the previous cases at 450 K. The 
H2 criterion for star formation is reached faster because 
of the excess electrons, which creates a star particle at a 
lower density (360 cm~^). 

4. Star formation after reincorporation (Halo 3, 4, 5) — 
After a sufhcient amount of gas that was expelled by dy- 
namical feedback of the first star is reincorporated into 
the halo, star formation is initiated again. Here virial 
temperatures of the halos are under 10^ K, which are 
hosting their second instance of star formation. These 
halos have a larger spread in gas densities and temper- 
atures than the halos forming their first star. Gas is 
heated by virialization and prior stellar radiation to over 
10^ K outside 10 pc. The central densities in halos 3 
and 5 are similar to the regions described in the first star 
formation section, however they are slightly warmer at 
410 and 480 K. Halo 4 shows a more diffuse core with 
densities of 550 cm~^ and temperatures of 440 K. 

5. Star formation in W^K halos (Halo 10, 11, 12) — 
In these halos, H2 formation is aided by atomic hydro- 
gen cooling. The ISM becomes increasingly complex as 
more stars form in the halo. The temperatures range 
from 100 K to 20,000 K throughout the halo. Halos 10 
and 12 have hosted 16 and 8 massive stars, respectively, 
since it started to continually form stars at z ~ 20. In 
SimB-RT (halo 10), the densities are higher than the 
cases. The gas in this halo is more centrally concen- 
trated than the others because the H II regions did not 
breakout of the halo, thus minimizing any outflows from 
feedback and dispersion of the central core. The tem- 
perature in halo 10 is significantly warmer than other 
regions at 2500 K. In halo 11, the devastation caused 
by three stars and their SNe at z = 21 prevented star 
formation until z = 16.9. Its initial recovery from that 
event is apparent by the single cool core with a tempera- 
ture of 590 K that sharply transitions to a warm, diffuse 
medium at r 10 pc. Halo 12 (SimA-RT) has a com- 
plex morphology that is not centrally concentrated and 
is caused by stellar outflows during a major merger (see 
iWise fc AbelllMFra for images). This morphology mani- 
fests itself in the radial profiles as large density contrasts 
spanning nearly six order of magnitude at r = 30 — 300 
pc. Similarly, temperatures range from 50 K to 10,000 
K in the same region. The star forms in a diffuse region 
(p = 130 cm^'^) that has a temperature of 470 K and 
whose H2 formation is enhanced because it resides in a 
relic H II region. 

4. STARTING COSMOLOGICAL REIONIZATION 

In this section, we first describe the ionizing radiation 
from massive stars that start cosmological reionization 
in small overdense regions we simulate. Then we discuss 
the effects of recombinations in the inhomogeneous IGM 
and kinetic energy feedback from Pop III stars. Lastly 
the evolution of the average IGM thermal energy is ex- 
amined. 

To illustratively demonstrate radiative feedback from 
massive stars on the host halos and IGM, we show pro- 
jections of gas density and temperature that are density- 
squared weighted in Figures [6] and [7] for all of the sim- 
ulations at redshift 17. These projections have the 
same field of view of 8.5 proper kpc and the same color 



maps. The large-scale density structure is largely un- 
changed by the stellar feedback, and the adjacent fila- 
ments remain cool since they are optically thick to the 
incident radiation. H2 cooling produces more centrally 
concentrated objects; however stellar feedback photo- 
evaporates <, IO^A/q halos near other star-forming ha- 
los. This is apparent in the density projections in 
the Jeans smoothin g around the most ma ssive halo (cf. 
iHaiman etlHI l2001l : iMesigner et all l2006f ). Kinematic 
feedback from SNe has an even larger effect on the sur- 
rounding gas. In SimB-SN, this effect is seen in the re- 
duced small-scale structure and low-mass halos with no 
gas counterparts. However, the most apparent difference 
in the radiative simulations is the IGM heating by Pop 
HI stars, especially in SimB-SN. 

4.1. UV Emissivity 

A key quantity in reionization models is volume- 
averaged emissivity of ionizing radiation. We utilize 
the comoving SFR to calculate the proper volume- 
averaged UV emissivity 



in units of ionizing photons per baryon per Hubble time. 
Here Qhi is the number of ionizing photons emitted in 
the lifetime of a star per solar mass, ~ 2 x 10^^ is the 
comoving mean number density, and 




30 25 20 15 



Redshift 

Fig. 8. — (a) Averaged emissivity in units of ionizing photons per 
baryon per Hubble time that is calculated from the star formation 
rate in Figure|3] (b) Mass-averaged ionization fraction (/e > 10"'^) 
of the inner 250 (300) comoving kpc for SimA (SimB). (c) Volume- 
averaged ionization fraction for the same runs, (d) The ratio of the 
mass- and volume-averaged ionization fraction. 
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stellar masses greater than 100 Mq, Qhi ~ 10 pho- 
tons per solar mass, corre sponding to 840 00 ionizing pho- 
tons per stellar proton (jSchaereil |2002| ). We plot the 
emissivity e in Figure [8^. It follows the same behav- 
ior as the SFR, but now can be directly used in semi- 
analytic reionization models. The emissivity increases 
from unity at redshift 30 to '--^lOO at the end of our simu- 
lations. Our results agree with the emissivity calculated 
in semi-analytic models that in clude Pop III stars (e.g. 
lOnken fc Miralda-Escudil2004[ ) and should be an upper 
limit however. 

4.2. Effective Number of Ionizations per UV Photon 

We show the mass-averaged and volume-averaged ion- 
ization fraction /e within the refined region in Figures[8)3 
and [SJ:. The ratio of these fractions are plotted in Fig- 
ure [SJi. The first star in the simulation ionizes between 
5-10% of the volume where we allow star formation. As 
stars begin to form in other halos after redshift 25, the 
ionization fraction gradually builds to 30% in the RT 
simulations and 75% in the SN case. The higher stellar 
luminosities in SimB-SN, which can be seen in Figure[8^, 
and the additional outflows generated by SN blast waves 
cause this difference in /e. The ratio {fe)m/ {fe)v illus- 
trates that the ionized regions are overdense by a factor 
of ~1.25 at early times in the simulation, but then de- 
creases to unity as the H II region grows. Additionally, 
the H II regions in halos with sustained star formation in 
the RT simulati o ns do not fully breakout into the IGM. 
iKitavama et al.l (|2004D provided a useful approximation 
of the critical halo mass 

in which an ionization front (I-front) cannot escape. 
This approximation is valid for stellar masses between 
80 and 500 M©, redshifts between 10 and 30, and sin- 
gular isothermal spheres. Our simulations exhibit this 
same trait in which I-fronts only partially breakout from 
the host halo above this mass scale. 

This is not the case with SNe because previous SN 
blast waves can more effectively evacuate the surround- 
ing medium, thus increasing the chances of radiation es- 
caping into the IGM from later stars in the same halo. At 
z = 21, there is an example of this occurring with three 
stars forming in succession in the most massive halo. Af- 
ter the first star goes SN, a diffuse and hot medium is 
left behind, but the blastwave has not completely dis- 
rupted two other nearby condensing clumps. The radia- 
tion from the second star now docs not have to ionize its 
host halo and has an escape fraction of near unity. The 
same happens for the third star in this halo. This episode 
further ionizing SimB-SN from 40% to 60%. As a note 
of caution, these ionized fractions should not be consid- 
ered as cosmological a verage because th ey only sample 
a highly biased region. Illiev et aTl ([2006) showed that a 
simulation box size of ^--^30 Mpc is needed to make global 
predictions. 

To examine the strength of recombinations, we com- 
pare the total number of electrons in the volume to the 
total number of ionizing photons emitted in Figure [9l 
The ratio of these two quantities is the number of UV 
photons needed for one effective ionization initially. This 




Redshift 

Fig. 9. — Top panels: Total number of ionizing photons emitted 
(thick lines) and total number of electrons (thin lines) for simula- 
tions with cooling only (dotted), star formation only (solid), and 
supernovae (dashed) in the inner 250 and 300 comoving kpc for 
SimA (left) and SimB (right). The H II regions are completely 
contained in these volumes. Bottom panels: The ratio of total 
number of electrons to the total number of ionized photons emit- 
ted. 
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Fig. 10.— Clumping factor C = {p^) / {p)"^ for SimA (left) and 
SimB (right), comparing the cases of the adiabatic equation of state 
(dot-dashed), atomic hydrogen and helium cooling (dotted), star 
formation only (solid), and supernovae (dashed). The clumping 
factors for all diffuse (<5 < 100) gas are plotted as thin lines, and 
the clumping factors for ionized, diffuse gas, Cion, are plotted as 
thick lines. 

ratio is approximately 3/5 (1/3) after the first star dies. 
The values in simulation A are higher due to its smaller 
volume. This ratio then steadily decreases from recombi- 
nations in the relic H II region to a few percent when the 
next star forms. As stars begin to form regularly in the 
simulation, there are 4 (6) photons per sustained ioniza- 
tion. However this ratio drops by a factor of 5 in the RT 
simulations after z ~ 20 when the H II regions become 
trapped in the halo, thus reducing the available photons 
for ionizing the IGM. The effects of SNe as previously 
discussed maintains the ratio of 6 photons required per 
ionization as the most massive halo grows in mass. 

4.3. Clumping Factor Evolution 

Volume averaged recombination rates in an inhomoge- 
neous IGM scale with the clumping factor C — {p^) /{p)^, 
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where the angled brackets denote volume averaged quan- 
tities. The recombination rate for hydrogen, e.g., is sim- 
ply 

=GWeP.(l+^f, (4) 

where fcrec is the case B recombination rate for hydrogen 



at r « 10* K, and fe 



s/n is the ionization fraction. 



Both the increased recombinations in overdense regions 
and photon escape fractions lower than unity result in 
the high number of UV photons needed for one effective 
ionization that we see in our simulations. 

Figure [TO] compares the clumping factor in the adia- 
batic, cooling only, star formation, and supernovae cal- 
culations, which we separate into C for all gas (thin lines) 
and ionized gas (thick lines). Since we resolve the local 
Jeans length by at least 4 cells in all simulations, the 
clumping factor is not underestimated, given our assump- 
tions about gas cooling in each model. Recall that our 
simulations use the Zel'dovich approximation for the ini- 
tial conditions, which may lead to underestimating the 
clumping factor. We do not include dense {S > 100) 
gas, residing in filaments and halos, in the clumping fac- 
tor calculation as we restrict this analysis to the IGM 
because these overdensities are self-shielded from inci- 
dent ionizing radiation (e.g. IShapiro et al.l I1997L 120041 : 
llliev et al.ll200a iWhalen et al.l l2008'). 

The clumping factor in the adiabatic and cooling-only 
cases smoothly increases to ^--^5 at z = 17 from unity 
at z > 30. The clumping factor in the star forma- 
tion and supernovae simulations are smaller by 10-25% 
than the other simulations because the overdensities in 
the IGM are photoevaporated by the nearby stars. SN 
explosions disperse gas more effectively than radiative 
feedback alone in larger halos and can have a bigger im- 
pact on the clumping factor. At redshift 20, the three 
stars and their SNe energy in the most massive halo de- 
stroy the surrounding baryonic structures and reduce the 
clumping back to the values seen in non-radiative cases. 

We plot the clumping factor Cion in ionized regions 
above > 10^^ as the thick lines in Figure [TOl This 
value is most relevant for recombination rates. Around 
the first star in the simulation, the IGM is overdense 
and contains a larger amounts of clumps than the rest of 
the simulation. This increased dumpiness competes with 
photoevaporation caused by the nearby star to create 
Cion values that are comparable to C for all gas at z > 
25. Afterwards Cion is always smaller than C with a 
maximum decrease of ~25%. 

4.4. Kinetic Energy Feedback 

SN explosion energy and kinetic energy generated 
in D-type I-front play a key role in star formation 
in low-mass halos, which are eas ily affected due to 
their shallow potential w eU (e.g . iDekel fc 8113 119861: 



Hachnch 1995; Bromm fcToebT l2003t IWhalen et all 



2004HKitavama et al.ll2004l: iKitavama k Yoshidall2005D . 



The kinetic energy created by SNe is sufficient to expel 
the gas from these low-mass halos. For example, the 
binding energy of a lO^M© halo is only 2.8 x 10^" erg 
at z — 20, which is two orders of magnitude smaller 
than a typical ene rgy output of a pair instability SN 
(|Heger fc WooslevI [2002 ). For a T^w > 10^ K halo at 
the same redshift, it is 9.4 x 10^^ erg. With our chosen 



stellar mass of 170 Mq, it takes 3-4 SNe to overcome 
this potential energy. 

The shock wave created by the D-type I-front travels at 
a velocity Vg — 25 - 35 km for density gra dients (i.e. 
p(r) P C r~") with slopes b e tween 1.5 and 2.25 (iShu et al.l 
[200a IWhalen et all [2001 iKitevama et all [20041) . This 
velocity is the escape velocity for halos with masses 
greater than 3 x 10® M© at z = 15, which is an order 
of magnitude greater than the most massive halos stud- 
ied here. However less massive halos can contain these 
I-fronts because pressure forces slow the I-front after the 
star dies. 

Using the position of the shock wave when the star 
dies (eq. [5]) and energy arguments, we can estimate the 
critical halo mass where the material in the D-type I- 
front can escape from the halo by comparing the binding 
energy Ei, of the halo and kinetic energy in the shell. 
For most massive stars, the shock wave never reaches 
the final Stromgren radius. 
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before the star dies. Here A^hi is the ionizing photon rate 
of the star, and n/ is the average number density of gas 
contained in this radius. After the lifetime of the star, 
the shock reaches a radius 



R^ = 83( -) 

VsOkms-i/ 



2.7Myry' 



(6) 



where is the stellar lifetime (see also [Kitavama et al.l 
[200l . We can neglect isolated, lower mass (M < 3OM0) 
Pop HI stars whose shock wave reaches i?str within its 
lifetime. In this case, the I-front stops at i?str, and the 
shock wave becomes a pressure wave that has no associ- 
ated density contrast in the neutral medium ([Shu[[199^ . 
Thus we can safely ignore these stars in this estimate. 

Assume that the source is embedded in a single isother- 
mal sphere. The mass contained in the shell is 



VsP^ 



(7) 



that is the mass enclosed in the radius Rg in an isother- 
mal sphere, corrected for the warm, ionized medium 
behind the I-front. Here Vs is the volume contained 
in a sphere of radius Rs, and pi is the gas density 
of the ionized medium, whose typical number density 
is 1 cm^'^ f or stellar feedback from a massive pri- 
mordial st ar (IWhalen et al.[ [2004 [Kit avama et al.|[200l 
f^shida et~aL| l2007al: [Abel et al.l [2007). For massive 
stars {Mi, >, 3OM0), the mass of the central homogeneous 
medium is small (i.e. 10%) compared to the shell. We 
compensate for this interior mass by introducing the frac- 
tion ry, so the shell mass is simply 



{VLb/nM) Mvir rs 



(8) 



For these outflows to escape from the halo, the kinetic 
energy contained in the shell must be larger than the 
binding energy, which is 
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Fig. 11. — Maximum halo mass in which a D-typc ionization 
front can create outflows as a function of primordial stellar mass 
for shock velocities Vs of 20, 25, 30, and 35 km s~^. Here the 
fraction r) of mass contained in the shell is 0.9. 



Using equations ([6]) and ([8|) in this condition, we obtain 
the maximum mass 



M 



Mmax - 3.20 X 10' 



VlOOpc J V30 km s-i 



0.17 
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of a halo where the material in the shock wave becomes 
unbound, expelling the majority of the gas from the halo. 

In Figure [TTl we use the stella r lifetimes and ionizing 
luminosities from lSchaered ()2002f ) to calculate the critical 
halo mass for outflows for stellar masses 5 - SOOAf© and 
for shock velocities of 20, 25, 30, and 35 km s~^ with rj = 
0.9. For stellar masses smaller than 3OM0, the D-type I- 
front reaches the final Stromgren sphere and cannot expel 
any material from the host. Hence they are not plotted 
in this figure. For the more massive stars, the star dies 
before the D-type I-front can reach the Stromgren radius, 
thus being limited by i*. This maximum halo mass is in 
good agreement with our simulations as we see halos with 
masses greater than 5 x IO^Mq retaining most of their 
gas in the star formation only cases. However in larger 
halos, stellar sources still generate champagne flows, but 
this material is still bound to the halo and returns in tens 
of million years. 

4.5. Thermal Energy 

Thermal feedback is yet another mechanism how Pop 
HI stars leave their imprint on the universe. The ini- 
tial heating of the IGM will continue and intensif y from 
hig her S FRs at lower redshifts (e .g. Hcrnquist & S pringell 
[20031: iOnken fc Miralda-Escudlr2004j) . It is possible to 
constrain the reionization history by comparing tempera- 
tures in the Lya fores t to different reionization scenarios 
()Hui fc Haimanir2003( ). Temperatures in the Ly a forest 
are approximately 20,000 K at 2 = 3 — 5 (Schav eet al.l 
I2OOOI: IZaldarri aga et al.ll2001h . Although our focus was 
not on redshifts below 15 due to the uncertainty of the 
transition to the first low-mass metal-enriched (Pop II) 



stars, we can utilize the thermal data in our radiation 
hydrodynamical simulations to infer the thermal history 
of the IGM at lower redshifts. 

The excess energy from hydrogen ionizing photons over 
13.6 eV photo-heat the gas in the H II region. The mean 
temperature within H II regions in our calculations is 
~30,000 K. When the short hfetime of a Pop HI star 
is over, the H II region cools mainly through Comp- 
ton cooling off the cosmic microwave background. The 
same framework applies to SNe remnants as well. The 
timescale for Compton cooling is 



tc = 1.4 X 10^ 



l + z 
20 



(11) 



This process continues until the gas recombines, and 
Compton cooling is no longer efficient because of its 
dependence on electron fraction. Radiation preferen- 
tially propagates into the voids and leaves the adjacent 
filaments and its embedded halos virtually untouched. 
Hence we can restrict the importance of Compton cool- 
ing to the diffuse IGM since Compton cooling cools the 
gas to low temperatures without being impeded by re- 
combinations that are proportional to n^. This causes 
the relic H II region to cool to temperatures down to 
300 K. The temperature evolution in our radiative cal- 
culations agrees with the analytic models of relic H II 
(|0h & Haiman 2003). 

We plot the volume- averaged temperature (T)v and 
mass- averaged temperature (T)ni in the volume where 
we allow star formation, i.e. the inner 250 (300) comov- 
ing kpc, in Figures [T^ and [T51 We compute the average 
temperatures in both neutral and ionized (/e > 10"'^) re- 
gions. We first focus on the thermal evolution of neutral 
gas. An increase of the average temperature in neutral 
gas is indicative of heating by hard photons or super- 
novae. At the final redshift z = 16.8, (T)v = 180 K 
in neutral gas with SNe compared with 90 K without 
SNe. Both of these average temperatures are a factor 
of 2-3 higher than without star formation. With SNe, 
(r)m also increases by -100 K to 1000 K. In the SimA 
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Fig. 12. — Evolution of the volume-averaged temperature in 
the inner 250 and 300 comoving kpc for Sim A (left) and Sim B 
(right), respectively, for neutral (thin) and ionized (thick; fe > 
10""^) gas. The simulations for the adiabatic (dot-dashed), cooling 
only (dotted) , star formation only (solid) , and supernovae (dashed) 
simulations are plotted. 
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panel of Figure I13[ radiative cooling in neutral gas by 
Lya and H2 in SimA-HHe and SimA-RT, respectively, 
is apparent at the final redshift. In the star formation 
run, most of the neutral mass fraction lie within the most 
massive halo and are shielded from radiation from stars 
within the halo and thus can radiatively cool. This oc- 
curs in SimB but to a lesser extent and is not clearly seen 
in Figure [T3l 

The effects from radiative feedback is most evident in 
the average temperatures of ionized gas. The first star 
in the calculations raises the mass- and volume-averaged 
temperatures of the ionized gas to 2 x 10"* K. Afterwards 
the remnant cools from Compton and adiabatic processes 
as it expands to temperatures similar to the RT simula- 
tions. Photoheating from later stars cause the temper- 
atures in the ionized regions to fiuctuate between 3000 
and 10000 K. The supernovae calculations are slightly 
higher due to the hot SN bubble that has an initial tem- 
perature of '^lO* K. The mass-averaged temperature in- 
creases more than (T)v because of the photo- heating of 
the host halo and virial heating of the halos, which is the 
cause of the increase in the simulations without star for- 
mation. These increased temperatures cause the photo- 
evaporation and Jeans smoothing of the gas in the relic 
H II regions. We discuss these effects in the next section. 

5. DISCUSSION 

We have studied the details of massive metal-free star 
formation and its role in the start of cosmological reion- 
ization. We have treated star formation and radiation in 
a self-consistent manner, allowing for an accurate inves- 
tigation of the evolution of cosmic structure under the 
infiuence of early Pop III stars. Stellar radiation from 
these stars provides thermal, dynamical, and ionizing 
feedback to the host halos and IGM. Although Pop III 
stars are not thought to provide the majority of ionizing 
photons needed for cosmological reionization, they play a 
key role in the early universe because early galaxies that 
form in these relic H II regions are significantly affected 
by Pop III feedback. Hence it is important to consider 
primordial stellar feedback while studying early galaxy 
formation. In this section, we compare our results to 
previous numerical simulations and semi-analytic mod- 
els of reionization and then discuss any potential caveats 
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Fig. 13. — Same as Figure [121 but for the mass- averaged temper- 
ature. 
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Fig. 14. — The Jeans mass Mj and filtering mass Mp that can 
form bound objects. The squares denote the total mass of star 
forming halos in all three simulations. 



of our methods and possible future directions of this line 
of research. 

5.1. Comparison to Previous Models 
5.1.1. Filtering Mass 

One source of negative feedback is the suppression of 
gas accretion into potential wells when the IGM is pre- 
heated. The lower limit of the mass of a star forming 
halo is the Jeans filtering mass 
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where a and Mj are the scale fac tor and time depen- 
dent Je ans mass in the H II region (iGnedin fc Huilfl998l : 
iGnedin 2000b) . Additionally, the virial shocks are weak- 
ened if the accreting gas is preheated and will reduce 
the collisional ionization in halos with Tvir ^ 10'* K. To 
illustrate the effect of Jeans smoothing, we take the large 
H II region of SimB-SN because it has the largest ionized 
filling fraction, which is constantly being heated after 
z = 21. Temperatures in this region fluctuates between 
1,000 K and 30,000 K, depending on the proximity of the 
currently living stars. In Figure [TH we show the result- 
ing filtering mass of regions with an ionization fraction 
greater than 10^'^ along with the total mass of star form- 
i ng halos. 

IGnedinI (|2000bfl found the minimum mass of a star 
forming halo is better described by Mp instead of Mj. 
Our simulations are in excellent agreement for halos that 
are experiencing star formation after reincorporation of 
their previously expelled gas. The filtering mass is the 
appropriate choice for a minimum mass in this case as the 
halo forms from preheated gas. However for halos that 
have already assembled before they become embedded 
in a relic H II region, the appropriate minimum mass 
-^min is one that is re gulated by th e LW background 
(jMachacek et al .l 120011; IWise & AbeH I20 05D and photo- 
evaporation (e .g. Efstathiou 1992; Barkan a fc Loeblll999t 
iHaiman et al][2001t iMesigner et al.ii2006} ). This is evi- 
dent in the multitude of star forming halos below Mp- 
With the exception of star formation induced by SN 
blast waves or I-fronts, this verifies the justification of 
using Afniin and Mp for Pop HI and galaxy formation. 
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respectively, as a criterion for star forming hafos in semi- 
analytic models. 

5.1.2. Star Formation Efficiency 

Semi-analytic models rely on a star formation efficiency 
fi,, which is the fraction of collapsed gas that forms stars, 
to calculate quantities such as emissivities, chemical en- 
richment, and IGM temperatures. Low-mass halos that 
form a central star have /* ~ 10""^ whose value originates 
from a single 100 Mp, star forming in a dark matter halo 
of mass 10^ Ma. (lAbel et al.l 120021: |Bromm et_al] 1200^ 
lYoshida et al.ll2006[ ). Pop 11 star forming halos are usu- 
ally calibrated with star formation efficiencies from local 
dwarf and high-redshift starburst gal axies and are usu- 
ally on the ord er of a few percent (e.g. lTavlor et al.lll999l : 
lGnedinl[2000ah . 

This leads to the question: how efficient is star forma- 
tion in these high-redshift halos while explicitly consid- 
ering feedback? This is especially important when halos 
start to form multiple massive stars and when metallic- 
ities are not sufficient to induce Pop II star formation. 
The critical metalli city for a trans i tion to Pop II is still 
unclear. Recently, iJappsen et al.l (|2007a!i ) showed that 
metal line cooling is dynamically unimportant in diffuse 
gas until metallicities of 10~^ Zq. On the other hand, 
dust that is produced in SNe can genera te efficient cool- 
ing d own in dense gas with 10~^ Zq (jSchneider et al.l 
[20061) . If the progenitors of the more massive halos did 
not result in a pair-instability SN, massive star forma- 
tion can continue until it becomes sufficiently enriched. 
Hence our simulations can probe the efficiency of this 
scenario of massive metal-free star formation. It has 
also been suggested that the cosmological conditions that 
lead to the collapse of a metal-poor molecular cloud 
{Z/Zq Ri 10~^'^) may be more important than some crit- 
ical metallicity in determ ining the initial mass function 
of a given stellar system (jJappsen et al.l[2007b( l. 

We calculate with the ratio of the sum of the stellar 
masses to the total gas mass of unique star-forming halos. 
For example at the final redshift of 15.9 in SimA-RT, the 
most massive halo and its progenitors had hosted 11 stars 
and the gas mass of this halo is 1.8x lO^M©, which results 
in /* = 6.1 X lO^** for this particular halo. Expanding this 
quantity to all star forming halos, /*/10^^ = 5.6, 6.7, 7.4 
for SimA-RT, SimB-RT, and SimB-SN, respectively. We 
note that our choice of M^, = UOMq in SimB-SN in- 
creases by 70%. Our efficiencies are smaller than the 
isolated Pop HI case because halos cannot form any stars 
once the first star expels the gas, and 40 - 75 million 
years must pass until star can form again when the gas 
is reincorporated into the halo. 

By regarding the feedback created by Pop HI stars 
and associated complexities during the assembly of these 
halos, the f^, values of ~6 x 10^^ that are explicitly de- 
termined from our radiation hydrodynamical simulations 
provide a more accurate estimate on the early star for- 
mation efficiencies. 

5.1.3. Intermittent & Anisotropic Sources 

Our treatment of star formation and feedback produces 
intermittent star formation, especially in low-mass ha- 
los. If one does not account for this, star formation rates 
might be overestimated in this phase of star formation. 



Kinetic energy feedback is the main cause of this behav- 
ior. As discussed in sections [3T2] and [44l shock waves cre- 
ated by D-type I-fronts and SN explosions expel most of 
the gas in halos with masses <, 10^ Mq. A period of qui- 
escence follows these instances of star formation. Then 
stars are able to form after enough material has accreted 
back into the halo. Only when the halo becomes massive 
enough to retain most of the outflows and cool efficiently 
through Lya and H2 radiative processes, star formation 
becomes more regular with successive stars forming. 

The central gas structures in the host halo are usu- 
ally anisotropic as it is acquiring material through accre- 
tion along filaments and mergers. At scales smaller than 
10 pc, the most optically thick regions produce shad- 
ows where the gas radially behind the dense clump is 
not photo-ionized or photo-heated by the source. This 
produces cometary and so-called elephant trunk struc- 
tures that are also seen in local star forming regions 
and have been discussed in detail since iPottaschI (jj958). 
At a larger distance, the surrounding cosmic structure 
is composed of intersecting or adjacent filaments and 
satellite halos that breaks spherical symmetry. The fil- 
aments and nearby halos are optically thick and re- 
main cool and thus the density structures are largely 
unchanged. The entropy of dense regions are not in- 
creased by stellar radiation and will feel little negative 
feedback from a n entropy floor that only exists in the 
ionized IGM (cf. lOh fc Haimanl [2003h . Ray-tracing al- 
lows for accurate tracking of I-fronts in this inhomoge- 
neous medium. Radiation propagates through the least 
optically thick path and generates champagne flows that 
have been studied exten sively in the context of present 
day star formation fe.g. [Franco et al.l 119901: IChurchweO 
[200llShu et"all[200a [Arthur fc H oarc 20Q3). In the con- 
text of massive primordial stars, these champagne flows 
spread into the voids and are impeded by the inflow- 
ing filaments. The resulting H II regions have "butter- 
flv" m orphologies ([Abel et al.lll999ll2007t [Alvarez et all 
I2006at iMellema et al.1 [20061: lYoshida et al.l l2007aD . We 
also point out that sources embedded in relic H II largely 
maintain or increase the ionization fraction. Here the al- 
ready low optical depth of the recently ionized medium 
(within a recombination time) allows the radiation to 
travel to greater distances than a halo embedded in a 
completely neutral IGM. The H II regions become in- 
creasingly anisotropic in higher mass halos. We show an 
example of the morphology of a H II region near the end 
of the star's lifetime in a dark matter halo with mass 
1.4 X IO'^Mq in Figure dni 

5.2. Potential Caveats and Future Directions 

Although we have simulated the first generations 
of stars with radiation hydrodynamic simulations, our 
methods have neglected some potentially important pro- 
cesses and made an assumption about the Pop III stellar 
masses. 

One clear shortcoming of our simulations is the small 
volume and limited statistics of the objects studied here. 
However, it was our intention to focus on the effects 
of Pop III star formation on cosmological reionization 
and on the formation of an early dwarf galaxy instead 
of global statistics. The star formation only simulations 
(SimA-RT and SimB-RT) converge to the similar aver- 
aged quantities, e.g. ionized fraction, temperatures, star 



14 



WISE & ABEL 





Fig. 15. — Density (left) and temperature (right) slices of an anisotropic H II region in the most massive halo of SimB-RT. The star has 
Uved for 2.5 Myr out of its 2.7 Myr Hfetime. The field of view is 900 proper parsecs. 



formation rates, at the final redshift. The evolution of 
these quantities differ because of the limited number of 
stars that form in the simulations, which then causes the 
evolution to depend on individual star formation times. 
This variance should be expected in the small volumes 
that we simulate and should not diminish the significance 
of our results. 

We have verified even in a 2.5-cr peak that Pop III stars 
cannot fully reionize the universe, which verified previous 
conclusions that low-luminosity galaxies provide the ma- 
jority of ionizing photons. Furthermore, it is beneficial 
to study Pop III stellar feedback because it regulates the 
nature of star formation in these galaxies that form from 
pre-heated material. Further radiation hydrodynamics 
simulations of primordial star and galaxy formation with 
larger volumes while still resolving the first star forming 
halos of mass ~3 x lO^Af0 will improve the statistics 
of early star formation, especially in more typical over- 
densities, i.e. 1-tT peaks, some of which could survive to 
become dwarf spheroidal galaxies at z = 0. 

In this work, we treated the LW radiation field as 
optically thin, but in reality, H2 produces a non-zero 
optical depth above co lumn densities of 10^^ cm^^ 
(jPraine fc Bertoldilll996f ). Conversely, Doppler shifts of 
the LW lines arising from large velocity anisotropics and 
gradient may render H2 self-shielding unimpo rtant up to 
colum n densities of 10^*^ — 10^^ cm^-^ (Glov er" Brandl 
I2OOII ). If self-shielding is important, it will lead to in- 
creased star formation in low-mass halos even when a 
nearby source is shining. Moreover, H2 production can 
also be cat alyzed ahead of I-fronts (iRicotti et al.|[200ll : 
lAhn fc Sha piro 200^. In these halos, LW radiation will 
be absorbed before it can dissociate the central H2 core. 
On the same topic, we neglect any type of soft UV 
or LW background that is created by sources that are 
cosmologically nearby (Az/z ~ 0.1). A soft UV back- 



ground either creates positi ve or nega tive feedback, de- 
pending on its strength ( Mesigner et al . 2006), and a LW 
background i ncreases the minimum halo mass of a star- 
forming halo (iMachacek et al.ll200 1':'Yoshida et al.ll2003l : 
lO'Shea fc NormanI 120071 : iWise fc Abel 2007b). However 
in our calculations, the lack of self-shielding, which sup- 
presses star formation in low-mass halos, and the neglect 
of a LW background, which allows star formation in these 
halos, may partially cancel each other. Hence one may 
expect no significant deviations in the SFRs and reion- 
ization history if one treats these processes explicitly. 

To address the incident radiation and the resulting UV 
background from more rare density fluctuations outside 
of our simulation volume, it will be useful to bridge the 
gap between the start of reionization on Mpc scales to 
larger scale (10 - 100 Mpc) simulations of reionization, 
such as the work of [Sp kasian ct al. (2003), Ilicv ct al] 
(jloni),[ZaMiEaD ([2007 '). and KohlereLaLI (|2007l) . Ra- 
diation characteristics from a volume that has similar 
overdensities as our Mpc-scale simulations can be sam- 
pled from such larger volumes to create a radiation back- 
ground that inflicts the structures in our Mpc scale sim- 
ulations. Inversely, perhaps the small-scale evolution of 
the clumping factor, filtering mass, and average temper- 
ature and ionization states can be used to create an ac- 
curate subgrid model in large volume reionization simu- 
lations. 

Another potential caveat is the continued use of pri- 
mordial gas chemistry in metal enriched regions in the 
SN runs. Our simulations with SNe give excellent ini- 
tial conditions to self-consistently treating the transition 
to low-mass star formation. In future work, we plan 
to introdu ce metal-line and dust cooling models (e.g . 
from Glo ver fc Jappsenl[2007l : iSmith fc SigurdssonllBoTl ) 
to study this transition. 

The one main assumption about Pop HI stars in our 
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calculations is the fixed, user-defined stellar mass. The 
initial mass function (IMF) of these stars is largely un- 
known, therefore we did not want to introduce an un- 
certainty by choosing a fiducial IMF. It is possible to 
calculate a rough estimate of the stellar mass by com- 
paring the accretion rates and Kelvi n-Helmholtz time 
of the contracting molecular cloud (|Abel et al.l 120021 : 
lO'Shea et al.l [20051 ). Protostellar models of primordial 
stars have also shown that the zero-age main sequence 
(ZAMS) is reached at 100 Mq for typical accretion 
histories after the star halts its adiabatic contraction 
(|Omukai fc Fallal 120031 : lYoshida et all l2006f l . Further- 
more, we have neglected HD cooling, which may be- 
come important in halos embedded in relic H II re- 
ions and r esult in lower mass {^30 Mp,) mctal-frc e stars 
O'Shea et al. 2005: Grcif & Bromm 2006; Yoshida,et_alJ 
120071)1) . Based on accretion histories of star forming ha- 
los, one can estimate the ZAMS stellar mass for each halo 
and create a more self-consistent and ab initio treatment 
of Pop III star formation and feedback. 

6. SUMMARY 

We conducted three radiation hydrodynamical, adap- 
tive mesh refinement simulations that supplement our 
previous cosmological simulations that focused on the 
hydrodynamics and cooling during early galaxy forma- 
tion. These new simulations concentrated on the for- 
mation and feedback of massive, metal-free stars. We 
used adaptive ray tracing to accurately track the re- 
sulting H II regions and followed the evolution of the 
photo-ionized and photo-heated IGM. We also explored 
on the details of early star formation in these simulations. 
Theories of early galaxy formation and reionization and 
large scale reionization simulations can benefit from the 
useful quantities and characteristics of the high redshift 
universe, such as SFR and IGM temperatures and ion- 
ization states, calculated in our simulations. The key 
results from this work are listed below. 

1. SFRs increase from 5x lO"'^ at redshift 30 to 6x 10"^ 
Mq yr-i Mpc~3 at redshift 20 in our simulations. Af- 
terwards the SFR begins to have a bursting nature in 
halos more massive than lO^M© and fluctuates around 
10^^ Mq yr^^ M pc~^. These rates are l arger than the 
ones calculated in iHernguist fc S pringel' ()2003[ ) because 
our simulation volume samples a highly biased region 
that contains a 2.5-cr density fluctuation. The associ- 
ated emissivity from these stars increase from 1 to ~100 
ionizing photons per baryon per Hubble time between 
redshifts 15 and 30. 

2. In order to provide a comparison to semi- analytic 
models, we calculate the star formation efficiency to be 
~6 X 10~^ averaged over all redshifts and the simula- 
tion volume. For Pop III star formation, this is a factor 
of tw o lower t han stars th at are no t affected b y feed- 
back ([Abel et al . 2002; Bro mm et al.ir2 002: YoshidaHTaD 
l2006t lO'Sheafc Normanll20'07ir 

3. Shock waves created by D-type I-fronts expel most 
of the gas in the host halos below ^5 x IO^A/q. Above 
this mass, significant outflows that are still bound to the 
halo are generated. This feedback creates a dynamical 
picture of early structure formation, where star forma- 
tion is suppressed in halos because of this baryon de- 
pletion, which is more effective than UV heating or the 



radiative dissociation of H2 . 

4. We see three instances of induced star formation in 
halos with masses ~ 3 x IO^Mq. Here a star forms as 
a SN blast wave overtakes an overdensity created by an 
ionization front instability. H2 formation is catalyzed by 
additional free elec trons in the r elic H II region and in 
the SN blast wave ()Ferraralll99"8[) . 

5. As star formation occurs regularly in the simulation 
after redshift 25, four (six) ionizing photons are needed 
per sustained hydrogen ionization. As the most mas- 
sive halo becomes larger than ^lO^Mg in the simulations 
without SNe, H II regions become trapped and ionizing 
radiation only escapes into the IGM in small solid angles. 
Hence the number of photons per effective ionization in- 
creases to 15 (50). In SimB-SN, stellar radiation from 
induced star formation have an escape fraction of nearly 
unity, which occur four times in the calculation. This 
allows the IGM to remain ionized at a volume fraction 
3 times higher than without SNe. Similarly, the ioniz- 
ing photon to ionization ratio also stays elevated at 10:1 
instead of decreasing in the calculations with star forma- 
tion only. 

6. Our simulations that include star formation and 
H2 formation capture the entire evolution of the clump- 
ing factor that is used in semi-analytic models to calcu- 
late the effective enhancement of recombinations in the 
IGM. We showed that clumping factors in the ionized 
medium fluctuate around the 75% of the values found in 
adiabatic simulations. They evolve from unity at high 
redshifts and steadily increase to ~4 and 3.5 with and 
without SNe at z = 17, respectively. Photo-evaporation 
from stellar feedback causes the decrease of the clumping 
factor. 

7. We calculated the Jeans filtering mass with the 
volume-averaged temperature only in fully and partially 
ionized regions, which yields a better estimate than the 
temperature averaged over both ionized and neutral re- 
gions. The filtering mass depends on the thermal his- 
tory of the IGM, which mainly cools through Comp- 
ton cooling. It increases by two orders of magnitude to 
~3 X IO^'Mq ai z ^ 15. It describes the minimum mass a 
halo requires to collapse after hosting a Pop HI star. For 
halos forming their first star, the minimum halo mass is 
regulated by the LW backgr ound (iMachacek et al.ll2001l) 
and photo-evaporation (e.g. iHaiman et al.ll2001l ). 

Pop III stellar feedback plays a key role in early star 
formation and the beginning of cosmological reionization. 
The shallow potential wells of their host halos only am- 
plify their radiative feedback. Our understanding of the 
formation of the oldest galaxies and the characteristics 
of isolated dwarf galaxies may benefit from including 
the earliest stars and their feedback in galaxy forma- 
tion models. Although these massive stars only partially 
reionized the universe, their feedback on the IGM and 
galaxies is crucial to include since it affects the char- 
acteristics of low-mass galaxies that are thought to be 
primarily responsible for cosmological reionization. Har- 
nessing observational clues about reionization, observa- 
tions of local dwarf spheroidal galaxies, and numerical 
simulations that accurately handle star formation and 
feedback may provide great insight on the formation of 
the first galaxies, their properties, and how they com- 
pleted cosmological reionization. 
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